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Abstract 

We have developed a methodology to study the thermodynamics of order-disorder transformations 
in n-component substitutional alloys that combines nonequilibrium methods, which can efficiently 
compute free energies, with Monte Carlo simulations, in which configurational and vibrational 
degrees of freedom are simultaneously considered on an equal footing basis. Furthermore, by ap- 
propriately constraining the system, we were able to compute the contributions to the vibrational 
entropy due to bond proportion, atomic size mismatch, and bulk volume effects. We have applied 
this methodology to calculate configurational and vibrational contributions to the entropy of the 
NisAl alloy as functions of temperature. We found that the bond proportion effect reduces the 
vibrational entropy at the order-disorder transition, while the size mismatch and the bond pro- 
portion effects combined do not change the vibrational entropy at the transition. We also found 
that the volume increase at the order-disorder transition causes a vibrational entropy increase 
of 0.08 kB/atom, which is significant when compared to the configurational entropy increase of 
0.27 kB/atom. Our calculations indicate that the inclusion of vibrations reduces in about 30% the 
order-disorder transition temperature determined solely considering the configurational degrees of 
freedom. 

PACS numbers: 63.50. Gh,81.30.Hd,65.40.gd,61.43.Bn 
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Monte Carlo simulations 



1 



I. INTRODUCTION 



One of the goals of materials science in the field of alloys is to predict and understand 
the relative stability of phases characterized by different chemical disorder. The disorder in 
an alloy can be considered as having a configurational contribution (configurational degrees 
of freedom), which is the disorder associated to the way the atoms are distributed in the 
parent lattice; and a vibrational contribution (vibrational degrees of freedom), which is the 
disorder associated to the phase space region around a static lattice configuration. For a 
very long time, most theoretical phase diagram calculations were done considering only the 
configurational degrees of freedorrti"^. In the 1990's, however, several experiments measuring 
thermodynamical properties of alloys in disordered metastable states,-'^'^'^'^'^'^ demonstrated 
the existence of a strong interplay between vibrational and configurational degrees of free- 
dom. It became clear that neglecting vibrational contributions to the thermodynamical 
properties of alloys could lead to inaccuracies, such as differences up to 30% between order- 
disorder (OD) transition temperatures calculated with and without vibrational degrees of 
freedom.— Theoretical studies of these alloys in a metastable disordered phase were per- 
formed assuming the alloys to be either completely ordered or totally disordered.— '^'^'^'^ 
The cluster variation method^ and its extensions^ have been used to calculate the config- 
urational contribution in partially disordered systems in equilibrium. Different approaches 
have been used to incorporate the vibrational degrees of freedom when cluster expansions 
are used, such as molecular dynamics,J^ the coarse graining method,— and the structure- 
inversion approach.— In the last two methods, the vibration contribution is taken into 
account through the harmonic approximation and anharmonicities are included via the 
quasiharmonic approximation. These two cluster expansion methods allow a first-principles 
description of the system, however, even calculations within the harmonic approximation 
are still very demanding for today's computer capabilities, and approximate approaches are 
still very useful.— These recent calculations^ have shown that anharmonic effects play an 
important role in the vibrational contribution to the thermodynamics of alloys. About 10 
years ago, a methodology that became known as MCX was proposed to treat simultaneously 
configurational and vibrational degrees of freedom by means of Monte Carlo simulations, 
which allow both atomic interchanges and atomic displacements^'^^'^i. 

In this work, we present a methodology to investigate phase equilibria of alloys that takes 
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into account naturally and simultaneously all configurational and vibrational contributions, 
including all anharmonicities, through a combination of the MCX approach and efficient 
tools to determine free energies, namely, the adiabatic switching (AS)^ and the reversible 
scaling (RS)^ methods. An interesting feature of our methodology is that it allows to study 
the contributions from different vibrational mechanisms to the total vibrational entropy. 
In real experiments, it is impossible to isolate the many factors that contribute to the vi- 
brational entropy such as the bond proportion, the atomic size mismatch, and the bulk 
volume mechanisms. On the other hand, it can be done in computer simulations, partic- 
ularly through the Monte Carlo technique, which makes possible to isolate the constraints 
associated with each mechanism. For example, the configurational contribution to the en- 
tropy can be calculated by allowing only the atomic interchanging dynamics. The effect 
of the bonds between different atomic species can be simulated by constraining the atoms 
to vibrate around their ideal crystalline structure positions and allowing the interchanging 
dynamics. The effect of the atomic size mismatch and the bond proportion mechanisms 
combined can be simulated by letting the atoms vibrate around their relaxed equilibrium 
positions at fixed volume and allowing the configurational dynamics. Finally, the volume 
mechanism can be simulated by allowing the volume of the supercell to vary by imposing 
constant pressure on the system, in addition to the positional and configurational dynamics. 
In summary, our methodology allows also to assess the contribution of a given vibrational 
mechanism by setting the appropriate constraint in the dynamics and then calculating the 
vibrational entropy difference between the relaxed and unrelaxed system. We have applied 
this methodology to quantify the vibrational entropy difference at the thermodynamical OD 
transition of the NisAl binary alloy. 

We chose the technologically important^i^^ NisAl as the alloy model for our study 
mainly because it is supposed to have the largest vibrational entropy difference upon 
disorder.-'^'^'^'^'^'^ The vibrational entropy difference due to disorder at the thermo- 
dynamical OD transition should be large enough to be unambiguously detected, since it 
is, in general, a fraction of the corresponding configurational entropy difference, which is 
itself relatively small. We also chose NisAl because it is particularly suitable to assess the 
magnitude of the size mismatch effect, since the difference between the atomic volumes of 
Al and Ni is quite large^^ {Vai - Vm) /{Vai + Vm) /2 = 0.41. In the case of NigAl, most of the 
research in the field, both experimental^'^ and theoretical, either using the embedded atom 
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method^'^'^ or tight-binding^^-^ potentials, has found a significant vibrational entropy dif- 
ference between the totally disordered (metastable) and the ordered phases. However, the 
subject is not free of controversy, van de Walle et alr^^, using ab initio calculations, found 
that the totally disordered and the ordered phases have essentially the same vibrational 
entropy. Our results indicate an increase of 0.08 kB/atom in the vibrational entropy at the 
thermodynamical OD transition, which is significant when compared to the corresponding 
configurational entropy increase of 0.27 ke/atom. We have also found that the bond propor- 
tion mechanism diminishes the vibrational entropy at the thermodynamical OD transition, 
whereas the size mismatch mechanism does not change it. These results are consistent with 
the local entropy calculations of Morgan et al.—'^. Regarding the importance of the volume 
mechanism, theoretical studie^^'^'^'^'^ have found that the volume mechanism is the main 
responsible for the increase of the vibrational entropy difference between the totally disor- 
dered and the ordered phase. This is supported by experimental work,— "^"^ in which it has 
been found an increase of the volume as the system becomes totally disordered. We have 
found that the volume increases 1.2% at the OD thermodynamical transition. In addition, 
our calculations indicate that the volume mechanism is the responsible for the increase in 
the vibrational entropy difference at the OD transition. 

The paper is organized as follows. In section II we present the general methodology. 
Section III describes the details of the interatomic potential we have chosen to describe the 
Ni3Al alloy. In section IV, the methodology is applied to evaluate the configurational and 
vibrational entropy as function of the temperature and the contribution of each vibrational 
mechanism at the OD transition. In section V we summarize the results. 

II. METHODOLOGY 

A. Dynamics and vibrational mechanisms 

In real systems, the process of chemical disordering takes place mainly through the migra- 
tion of vacancies.—"^ The problem of realistically simulating the disordering process through 
this mechanism is that the average vacancy concentration is very low (less than 10~^),~ im- 
plying the requirement of very large system sizes. For this reason, we chose the atomic 
exchange dynamics to simulate the chemical disorder. This dynamics can be implemented 
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through the Monte Carlo method. In this approach, the configurational degrees of freedom 
are explored by selecting at random two atoms belonging to different chemical species, their 
positions in the lattice are then interchanged, the energy change upon the atomic exchange 
is calculated, the Boltzmann factor associated to this change in energy is computed, and 
the move is accepted or rejected according to the Metropolis algorithm.— In order for the 
algorithm to be efficient, one should keep two lists of atoms of each atomic species and 
choose randomly one atom of each list to form the pair of atoms to be interchanged. This 
can be easily done, since the number of atoms of chemical species is kept constant. This 
dynamics is more efficient than the Kawasaki dynamics and satisfies detailed balance. We 
will call this dynamics as the configurational case. We consider a Monte Carlo step (MCS) 
in the configurational case as attempts to exchange the atomic positions of two atoms of 
different species chosen at random, where is the number of atoms. 

In order to investigate the various contributions to the vibrational entropy we considered 
different dynamics related to the various mechanisms we have mentioned earlier. To esti- 
mate the effect of the bond proportion mechanism in the vibrational entropy we define the 
following dynamics: besides the atomic interchange dynamics, we constrain the atoms to 
vibrate around their ideal crystalline structure positions. We call this dynamics as the un- 
relaxed case. The vibrational dynamics for the unrelaxed case is accomplished by choosing 
one atom at random and calculating a new state by 

xr"' = xf + A^,,(26-1), t = 1..3, (1) 

where xf is the coordinate associated to the corresponding ideal crystalline structure posi- 
tion, is a random number between zero and one, and Amax is the maximum displacement 
allowed, which is adjusted automatically in such a way that 50% of the trials are accepted.— 
In this case, a MCS was considered as attempts of atomic displacements followed by N' 
attempts to exchange atoms. We chose A^' = A^/10, in the particular case of NisAl, because 
it is the minimum number of attempts of atomic exchange needed for the potential energy 
and the order parameter to relax to average values. The bond proportion mechanism is 
claimed^'^'^ to be relevant for changes in the vibrational entropy due to chemical disorder, 
since the proportion of bonds between distinct and similar atoms changes with the disor- 
der. This effect could in principle be addressed, from the theoretical point of view, through 
calculations based on a spring model,— "^^"^ which presumably assigns a softer and a stiffer 
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character to the bonds between like and unlike atoms, respectively.— 

We now set a different type of dynamics by letting the atoms relax around their thermal 
equilibrium position as 

This kind of dynamics allows the atoms to vibrate around their thermal equilibrium posi- 
tions, taking into account naturally their accommodation due to their different atomic sizes. 
This is another mechanism that is claimed to influence the vibrational entropy and is called 
the size mismatch mechanism. It can be understood by the following picture;^ when two 
large atoms are constrained in a small space they can experience a compressive stress, in- 
creasing the stiffness of the bond, and reducing the allowed region to vibrate; on the other 
hand, when two small atoms are constrained in the same space they can experience a tensile 
stress, increasing the softness of the bond, and increasing the allowed region to vibrate. The 
interplay between the compressive and tensile stresses will dictate the final vibrational ef- 
fect. We call this relaxation of the configurational and positional constraints as the partially 
relaxed case. With this dynamics we simulate the combined effect of bond proportion and 
size mismatch in the vibrational entropy. In this case also, a MCS was considered to be the 
same as in the unrelaxed case. 

Finally, we simulate the contribution from the volume mechanism by allowing the volume 
change in a constant pressure simulation. This effect is explained^ by the following reason- 
ing: as the average interatomic distances increase, the bonds become softer and the atoms 
have more space to visit, that causes an increase in the vibrational entropy. We call this 
relaxation of the configurational, positional, and volumetric constraints as the fully relaxed 
case. Through this dynamics we simulate the combined effect of the bond proportion, size 
mismatch, and volume mechanisms in the vibrational entropy. The volumetric dynamics is 
accomplished by rescaling the atomic positions according to the volume change calculated 
as 

L-^- = L"'' + ^LM-^). (3) 

where L is the length of the simulation box, and A^„^ is the maximum allowed change 
in that length, which is adjusted similarly to the unrelaxed and partially relaxed cases at 
typically 10 MCS. In the fully relaxed case, a MCS was defined as positional trials followed 
by N' = N/10 exchanging trials, and one volumetric trial.— 
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B. Free energy calculations 

The thermodynamical quantity that underhes all this work is the free energy, which is 
calculated through the AS^ and RS^ methods. The AS method allows one to calculate the 
free energy by computing the work done by adiabatically switching the Hamiltonian of the 
system of interest to the Hamiltonian of a reference system (or vice- versa) at a single given 
temperature. On the other hand, the RS method allows one to evaluate the free energy 
in a range of temperatures provided that it is known at a single given temperature. These 
methods are very efficient since they evaluate the free energy from only one simulation run, 
whose length is determined by the required accuracy. In contrast with other methods, such 
as the harmonic,— or the quasi-harmonio^'^'^ approximations, the AS and RS, take into 
account naturally all anharmonic effects, which are crucial for the calculation of vibrational 
entropy differences. 

The AS method is based on the well known thermodynamic integration method (TI).— 
In the TI method, the absolute free energy of a system of interest can be estimated by 
computing the work done to transform the Hamiltonian of a reference system, of which one 
knows the free energy, into that of the system of interest. This can be achieved by considering 
the artificial Hamiltonian, H{X) = XHsys + {1 — X)Href, where H^ys is the Hamiltonian of the 
system of interest, H^ef is the Hamiltonian of the reference system, and A is a dimensionless 
coupling parameter. By varying A from to 1, one can transform one Hamiltonian into the 
other one. The work performed to switch between the two systems is given by the integral 



If now A is considered to be a function of time, and its value continuously varied from to 1 
during the time of simulation ts, the free energy difference between the two systems is given 



where Usys is the potential energy of the system of interest, Uref is the potential energy of 
the reference system, Wirr and Wrev are the irreversible and reversible work, respectively, 
and Ediss is the energy dissipation. Time in Eq.(l5]) can be regarded as the actual time, as 
in a molecular dynamics simulation, or the fictitious time created by the successive steps 
in a Monte Carlo simulation. The potential energy difference between the reference system 
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and the system of interest appears in Eq.([5|), instead of the Hamiltonian difference, because 
we consider the kinetic degrees of freedom to be in equilibrium, and, therefore, the kinetic 
energy terms cancel each other. The energy dissipation is one source of error, characteristic 
of nonequilibrium dynamic processes, and can be estimated^ by performing the direct and 
inverse transformations between the two systems: 

Tj^ref^sys , rr^sys^ref 

= ^ T^^- 

In all AS and RS calculations we adopted this criterion to quantify the free energy error, 
which can be reduced by increasing the simulation time. Another source of error-^ are the 
statistical fluctuations of the quantities in the integrand of Eq.Q, which can be handled by 
simulating other trajectories and averaging the results. 

There are subtleties in the AS method we must be aware of in order to obtain the correct 
free energy. The external conditions, like temperature, and the parameters of the reference 
system must be set in a way that the coupled system, described by H{X), does not undergo 
a phase transition along the transformation path. (We took particular care about the choice 
of the reference temperature for the high temperature disordered phases in the unrelaxed, 
partially relaxed, and fully relaxed cases to avoid mechanical melting.) 

Now let us discuss briefly the RS method and its application.— In contrast to the AS, 
the RS method allows the calculation of the free energy in a range of temperatures. This 
can be accomplished by realizing that the free energy of the scaled system at a temperature 
To, whose potential energy is given by U scaled = XUgys (in the case of RS A is not restricted 
to the interval [0, 1]), is related to the free energy of our system of interest at a temperature 
T = To/A.— 1^ The free energy of the scaled system at a given value of A can be readily 
determined by computing the work performed to change A from 1 to A = Tq/T, as it is done 
in the AS method (provided that the free energy is known for A = 1). It is shown in Ref. 



26l that the free energy of a system at temperature T can be estimated from the irreversible 



work Wirrit) done to bring the system from Tq to T(t), as 

Fim) F{n) W^t) 3 T(t) 

where T{t) = To/A(t), and T(To) is the known free energy reference. The logarithmic term 
of Eq.([7l) corresponds to the contribution of the kinetic degrees of freedom and must be 
omitted when only the conflgurational changes are considered. The estimation of the energy 
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dissipation at a given temperature is calculated, as in the AS method, using Eq.([6]). In the 
case where the external pressure is set zero, as in the calculation of the free energies for the 
fully relaxed case, the Gibbs free energy formula reduces to Eq.Q. 

III. THE NigAl SYSTEM 
A. The choice of potential 

Some well known and often used interatomic potentials for modeling NisAl^i^ are not 
suitable to describe the configurational degrees of freedom of this alloy.-^ The reason for 
that is that this potential, in both parameterizations,— does not yield the LI2 phase 
as the ground state phase, giving rise to nonphysical thermodynamical phases at low 
temperatures.— In order to modeling appropriately the NisAl system we looked for a poten- 
tial which provides not only the correct ground state, but describes, at least qualitatively, the 
thermodynamics of the OD and vibrational phenomena. We chose the tight-binding Finnis- 
Sinclair— potential whose parameterization was obtained by Vitek et al.— Among the ther- 
mal features of this potential we may cite the linear thermal expansion coefficient (at 1050 
K) of 21.7 X lO^^K^^, which agrees well with the experimental value of 19 x lO^^K"^;^ the 
equilibrium lattice parameter (at 1000 K) of oq = 3.6096 A, which is in good agreement with 
the experimental value of oq = 3.6120 and the calculated mechanical melting tempera- 
ture T^'^'^^ = 1600 ± 25 K at the Lindemanns'(5 function value of 0.12. We have determined 
the thermodynamical melting temperature for this model of NisAl to be = 1328 ± 6 
K. The thermodynamical melting point of a substance is obtained by determining at which 
temperature the solid and liquid phases have the same free energy. It is important to point 
out that the thermodynamical melting temperature we have obtained for the model is 20% 
lower than the experimental value = 1636 K.— This discrepancy in the melting transi- 
tion temperature is not surprising since the potential parameters are fitted from a database 
which does not include data from the liquid phase. We will return to this point later, after 
we present the results for the OD transition. However, the important conclusion we should 
advance at this point is that, despite numerical discrepancies, the results from our simula- 
tions for the OD transition temperature and the melting point are qualitatively consistent 
with experimental findings. 
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B. Order parameters 



In the case of LI2 alloys, the order parameters can be defined as follows. The long-range 
order parameter is constructed from the LI2 phase by labeling the sublattice associated to 
the Al (Ni) atom as an a (/?) sublattice. The long-range order is then measured by the 
formula, first introduced by Bragg and Williams,- 

Pa - 0-25 

where Pa means the fraction of Al atoms in the a sublattice. In this way, rj = 1 for the 
ordered LI2 phase and r/ = for the disordered phase. This order parameter is very useful 
to quantify the long-range order, however one must be careful with its interpretation. First, 
when one performs computer simulations to explore the configurational degrees of freedom 
through cooling experiments, at a relative low rate, starting from the disordered phase at 
high temperatures, the system should always end up in the LI2 phase, but the Al atoms not 
always are found in the arbitrarily defined a sublattice. In other words, one does not know, 
in advance, which one of the four possible sublattices will be a sublattice. Hence, in this 
kind of experiment one must measure the long-range order parameter in the four possible 
sublattices. Second, when the system is in an antiphase boundary (APB)-like configuration^ 
a large fraction of Al atoms may be in a ordered block at sites of a /5 sublattice, giving low 
and even negative values for rj. Therefore, r] = does not distinguish between a totally 
disordered and a particular APE configuration. In Fig. [Dd, we show the results for the 
long-range order parameter for the a and /3 sublattices. 

Concerning the short-range correlations we measure the short-range order parameter, 
first introduced by Bethe and Wills,-^ as 

PAl-Ni - 9 



a 



12-9 ' 



(9) 



where pAi-Ni means the average number of unlike bonds between an Al atom and its first- 
neighbors. Note that a = implies r] = 0, however, t] = can correspond to cr = 1 as in a 
particular APE configuration. 
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C. Implementation details 



We have performed tests of our computational code by calculating the melting temper- 
ature of the Ni system using the Cleri and Rosato potential,— which agrees exactly with 
the value reported in Ref. l53|. Furthermore we compared our Finnis-Sinclair calculations of 
antiphase-boundary and stacking fault energies with Vitek et al.— results, and verified an 
exact agreement. 

The reference system chosen for the calculation of free energy references in the solid 
phases was the Einstein crystal.— "^"^ The chosen values for the vibration angular frequen- 
cies are wai = 75.4 rad THz, and w^i = 31.4 rad THz, for the Al atom and for the Ni 
atom respectively. These are the frequencies of the main vibrational modes of the elements, 



56l . which are expected 



estimated from the experimental phonon density of states from Ref. 
to be optimal to mimic the vibration of the atoms in the alloy. The reference system chosen 
for the calculation of the free energy reference of the liquid phase (used to estimate the 
melting temperature) was the r~^^ repulsive fluid.— The repulsive fluid parameters are 
chosen in such a way that the position and height of the flrst peak of the radial distribu- 
tion function of the r"^^ repulsive fluid potential coincide with those of the Finnis-Sinclair 
potential. This choice of parameters enhances the probability of the reference system to be 
within the borders of the phase diagram of the system of interest, thus minimizing the risk 
of encountering a phase transition.— "^"^ 

In the AS and RS calculations we chose time simulations such that the energy dissipation 
was less than 10""^ eV/atom, which typically leads to simulation lengths of 2 x 10^ MCS. 
The functional form of X(t) was always chosen to be a linear interpolation between the 
initial and flnal simulation times, which correspond to the initial and final temperatures. To 
circumvent surface effects we applied periodic boundary conditions and the minimal image 
convention.— Since both, the Einstein crystal and the r~^^ repulsive fluid do not have any 
cohesion, the simulations involving these systems have to be performed at flxed volume, 
which is chosen to be the average volume of a NPT equilibrium simulation at the given 
pressure and temperature of interest. Aside from the systematic errors due to dissipation, 
statistical errors in free energy calculations were handled by taking averages over typically 
10 samples. The error bars in the entropy differences were obtained from the fluctuation of 
entropy data below and above the transition in the standard way. In most of the results that 
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will be presented in the following section, a simple running average smoothing procedure was 
used in order to remove the unwanted fluctuations introduced by the numerical derivative 
calculations. All the calculations were performed using a cubic simulation cell containing 
500 atoms. Finally, we tested a larger system size using a 1372-atom simulation cell and 
found no signiflcant finite-size effects in entropy differences and transition temperatures. 



IV. RESULTS AND DISCUSSION 
A. The configurational case 

Let us discuss briefly the equilibrium numerical experiments performed in order to bracket 
the OD temperature for the configurational case. We set the LI2 phase at a fixed volume 
corresponding to the equilibrium volume obtained at zero pressure and Tq = 10^ K. Then 
we turned on the exchange dynamics and performed a series of equilibrium simulations on a 
relatively fine grid over a temperature interval of 2000 K. The measured thermodynamical 
quantities are shown in Fig. [H The abrupt changes in the behavior of the potential energy, 
specific heat, and order parameters indicate the OD transition at T^^"-^ = 1925 ±30 K. Note 
the abrupt change of the long-range order parameter 77, and the nonzero value of a after 
the transition. In order to estimate the effect of the chosen fixed volume on our results, 
we performed an analogous series of calculations at the equilibrium volume at K, which 
is substantially smaller than the one at zero pressure and Tq = 10^ K. We found the OD 
transition to be only 5% larger than the previous one, and such small difference indicates 
that the chosen value for the volume is not relevant for our conclusions. 

We now discuss the free energy calculations. We consider the free energy reference in 
this case to be at infinite temperature, because in this limit the configurational entropy per 
atom of a system with N atoms has an analytical expression corresponding to the ideal solid 
solution, which is 

1 

S^nfioo) = J^kslnj^-^. (10) 

This quantity measures the number of distinct configurations obtained by arranging N^i and 
N^i atoms in the lattice. Thus, for a system containing 500 atoms, we have approximately 
Sconfi^o) = 0.556 ke/atom. The advantage of using the RS method in this case results from 
the fact that the method maps the problem of determining the free energy for an infinite 
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interval of temperature onto a problem of finding the free energy for a finite interval of the 
scaling parameter A. We determine the work done to take the system from Tq = 10'^ K 
(A = 1) to the virtually infinite temperature (A = 0). Combining this work and the entropy 
from Eq. lfTOll . we are able to calculate the free energy at To using Eq.Q, noting that in the 
configurational case the logarithmic term should be dropped, since in this case the atoms 
are not allowed to vibrate. From the free energy reference at Tq and by computing the work 
done to take the system from A = 1 to any A < 1, we are able to calculate F as a function 
of T as shown in Fig. [2l 

In order to estimate the energy dissipation we compute the work performed to take the 
system from A = (infinite temperature) to A = 1 (Tq). As we can see in Fig. Eh the energy 
dissipation in the direct and reverse processes is less than 10""^ eV/atom. The configurational 
entropy, shown in Fig. [2b, is obtained by computing numerically —d{Fconf) / dT , where the 
brackets denote an average over uncorrelated samples of the configurational free energy. 
This averaging procedure is done in order to reduce the statistical noise in the numerical 
derivative (subsequently the remaining noise is further reduced by applying a running average 
smoothing procedure). The OD transition temperature, which in this case is considered to 
be at the center of the coexistence region of the two phases, was found to be T™"-^ = 1925 ± 
30 K, which agrees with the transition temperature determined by analyzing the behavior 
of thermodynamical quantities in Fig. [TJ The configurational entropy difference calculated 
at the OD transition is AS'^°^J = 0.18 ± 0.01 ke/atom. In Fig. [3b we show the behavior of 
the order parameters as function of the temperature. We can see that the long-range order 
parameter goes to zero at the OD transition temperature, whereas the short-range remains 
finite above the transition, approaching zero only for extremely high temperatures. Due to 
the persistence of the short-range order, we note that the configurational entropy reaches 
its maximum only at very high temperatures. The similar behavior of the configurational 
entropy and short-range order parameter with temperature allows us to establish a direct 
relationship between these two magnitudes. This relationship will be used in the calculation 
of free energy references for the other cases we will study. This result may also be especially 
useful at much higher pressures where the OD temperature is much lower than the melting 
point.— 
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B. Unrelaxed, partially relaxed, and fully relaxed cases 



Now we turn to the other cases, where the constraints on the degrees of freedom are 
gradually lifted. When vibrations are allowed, the limit of infinite temperature is no longer 
a suitable reference for the free energy, since the system would not remain a solid. In order to 
bracket the OD transition temperatures for each case, we performed a series of heating and 
cooling experiments. In Fig. [4] we plot the short-range order parameter a as a function of the 
temperature, for the unrelaxed, partially relaxed, and fully relaxed cases. The metastability 
exhibited in each case allows suitable choices of reference temperatures for the calculation of 
reference free energies using the AS method. For each case, two reference temperatures are 
chosen, one below and one above the guessed OD transition temperature provided by the 
metastability region. These free energy references are subsequently used in the RS method 
to calculate the free energy curves starting from both reference temperatures. Starting at 
the lower reference temperature, the RS method generates a free energy curve for increasing 
temperatures; from the higher reference temperature, on the other hand, the RS method 
provides a free energy curve for decreasing temperatures. The crossing of these two curves 
gives the OD transition temperature. These OD transition temperatures are indicated by 
dotted lines in Fig. [H The error bars for the transition temperatures were obtained from the 
free energy reference error bars in the same way as in Ref. |6ll: the RS method is performed 
again, this time starting from the extremes of the error bar of free energy reference given by 
the AS method. The temperature error bar is then obtained by the two crossings points that 
are the farthest from each other among the four curves intersections around the transition. 
Next, we discuss the details of these calculations and further results. 

The total free energy at the reference temperature T^g/ for the unrelaxed, partially re- 
laxed, and fully relaxed cases is calculated by adding the contributions from the vibrational 
and configurational degrees of freedom as 

F{Tref) = F^ih^Tref) — T^ef S conf{(^{Tref)) i (H) 

where FmbiTref) is the vibrational free energy calculated through the AS method to a refer- 
ence system which does not take into account the configurational entropy, e.g., the Einstein 
Crystal, and Sconf is the configurational entropy corresponding to the short-range order 
parameter at T^e/. The mapping between Sconf and cr is obtained from their temperature 
dependence in the configurational case, using the data given in Fig. [3l 
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In order to evaluate the contribution of the bond proportion mechanism to the vibrational 
entropy, we calculated the free energy of the system in the unrelaxed case. In Fig. [5] we 
show the free energy and entropy curves below and above the OD transition. The crossing 
of the two free energy curves determines the OD transition. The entropy is then computed 
by the numerical derivative of the free energy curves. The OD transition temperature 
obtained where the two curves cross is T^J = 1635 ± 60 K. Once the transition temperature 
is obtained, one can go back to the data in Fig. [4] to determine the short-range order 
parameter for the ordered and disordered phases at the OD transition, and from that the 
configurational entropy difference at the OD transition. The total entropy difference and 
the configurational entropy difference at the transition are AS**"* = 0.20 ± 0.02 ke/atom 
and AS!^°"'^ = 0.27 ± O.OlkB/atom, respectively. So the entropy difference due to only the 
vibration 

As:f = AS':: - A^r^ (12) 

is — 0.07±0.02kB/atom. This negative value is consistent with the local entropy calculations 
of Morgan et al.,—'^ who concluded that the increase of Al nearest-neighbors decreases the 
local vibrational entropy of Al or Ni atoms. 

The combined effect of size mismatch and the bond proportion mechanisms was studied 
by simulating the system in the partially relaxed case. The OD temperature obtained by the 
crossing of the free energy curves is T^J = 1497 ± 40 K. The total partially relaxed entropy 
difference and the configurational entropy difference at the transition are AS**"* = 0.23 ± 
0.02 ke/atom and AS^"""^ = 0.22±0.01 ke/atom, respectively. So the entropy difference due 
to only the relaxation of the atoms at fixed volume is approximately 0.01 ± 0.02 ke/atom. 
This result is also consistent with Morgan et al.^^ who observed that the local vibrational 
entropy "seems very fiat, or even slightly increasing" as the number of relaxed Al first- 
neighbors increases. 

Simulations of the system in the fully relaxed case were employed to evaluate the combined 
effect of the volumetric relaxation of the supercell with the size mismatch and the bond 
proportion mechanisms. The calculated OD temperature Tj^^ = 1339 ± 20 K essentially 
coincides, within the error bars, with the melting temperature of Tm = 1328 ± 6 K. This is 
in agreement with experimental findings.—"^ Although the Finnis-Sinclair potential cannot 
reproduce quantitatively the values obtained experimentally, it provides results that are 
consistent with the experimental results. 
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In order to show the significance of the configurational disorder on the vibrational prop- 
erties of the alloy, we depict in Figl6] the vibrational entropy of the alloy in the fully relaxed 
case as a function of temperature in comparison with the total entropy of the alloy in the 
fully relaxed case and in the perfectly ordered LI2 phase. In the fully relaxed case the 
vibrational entropy is obtained by subtracting the configurational entropy from the total 
entropy. The configurational entropy is, in turn, obtained as a function of temperature from 
its mapping with the measured short-range parameter. The entropy of the LI2 phase is 
purely vibrational. Therefore, the difference between the fully relaxed vibrational entropy 
and the LI2 vibrational entropy is only due to configurational disorder. We see that, as the 
temperature increases, the vibrational entropy due to disorder increases steadily and has a 
finite jump at the OD transition. 

The total entropy difference and the configurational entropy difference at the OD tran- 
sition in the fully relaxed case are ASjf = 0.35 ± 0.02 ke/atom and AS^J^-^ = 0.27 ± 
0.01 ke/atom, respectively. So the entropy difference excluding the configurational contri- 
bution is ASjf = 0.08 ± 0.02 ke/atom. This result shows that the vibrational entropy 
difference at the OD temperature is about 23% of the total entropy difference when the 
volume is allowed to relax. Furthermore, this vibrational entropy increase is accompanied 
by a volume increase of 1.2%. This result, together with the essentially zero vibrational 
entropy difference found in the partially relaxed case, indicates that the volume mechanism 
is the responsible for the vibrational entropy increase in NisAl. The increase in volume upon 
disorder is consistent with all experimental^'^'^'^ and theoretical work— The re- 
sult that the vibrational entropy difference increases at the OD transition is consistent with 
all the experimental^'- and most of the theoretical work,— '^'^'^'^ which have observed a 
positive vibrational entropy difference between the totally disordered (metastable) and the 
totally ordered phases. The OD temperature in the fully relaxed case is approximately 30% 
lower than the OD temperature calculated when only the configurational degrees of freedom 
are considered. Ozoli^s, Wolverton, and Zunger— proposed a relation between the OD tran- 
sition temperature calculated considering only the configurational degrees of freedom and 
the transition temperature determined including also vibrations. 



Using the results from our calculations as input for Eq. (fT3ll . namely, T^^"-' = 1925 K, 




(13) 
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AS'f^2j = 0-18 kB/atom, and AS]f> = 0.08 kB/atom, we find T™''/+«^ = 1333 K, i.e., the 
inclusion of vibrations lowers the OD transition temperature by 31%, with respect to the 
purely configurational transition temperature, which agrees with the lowering of 30% we 
have determined in the fully relaxed case. 

Finally, in order to compare the changes in entropy for all cases studied (and also for the 
liquid phase) we show in Fig. [7] the entropy as a function of the temperature. 

V. SUMMARY 

In this work we employ to the greatest possible advantage the RS and the Monte 
Carlo techniques, providing a methodology to evaluate both the configurational and vi- 
brational free energies as functions of temperature for n-component substitutional alloys. 
This methodology is also used to quantify the contributions of the vibrational mechanisms 
by setting appropriate constraints in the dynamics and calculating the entropy difference 
between the relaxed and unrelaxed system. By consecutively relaxing the configurational 
and structural constraints, we are able to quantify the configurational entropy as well as the 
vibrational entropy due to the bond proportion, size mismatch, and volume mechanisms. We 
applied this methodology to study the vibrational entropy difference at the thermodynamical 
OD transition of the NisAl alloy obtaining the following results. When the atoms are allowed 
to both interchange their positions and vibrate around their ideal fee lattice positions, the 
vibrational entropy difference is AS^f = — 0.07 ± 0.02 kB/atom. This indicates that bond 
proportion mechanism decreases the overall vibration of atoms upon transition. When the 
atoms can interchange their positions and are allowed to vibrate around their equilibrium 
positions, the vibrational entropy difference is essentially zero. This indicates that the size 
mismatch, coupled to the bond proportion mechanism, does not change the vibrational en- 
tropy upon transition. When the constraints on atomic interchanging, atomic relaxations, 
and bulk volume are relaxed, the vibrational entropy difference at the OD transition is 
ASjf = 0.08 lb 0.02 kB/atom, which is substantial when compared with the configurational 
entropy difference of ASj^J^^ = 0.27±0.01 kB/atom. This indicates that the effect of volume 
relaxation is the source of the increasing in the overall vibrational entropy upon disorder. 
Moreover, the volume increase at the OD transition is of 1.2%. A particular relevant result 
is that the OD transition temperature calculated when all constraints are allowed to relax 
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is approximately 30% less than that calculated when only the configurational degrees of 
freedom are considered. This result corroborates the importance of the vibrational degree 
of freedom in the determination of precise OD phase diagrams. Finally, as this methodology 
is not restricted to a particular crystal structure and stoichiometry, it can be applied to any 
n-component substitutional alloy to evaluate the configurational and vibrational entropies 
as function of the temperature and quantify the importance of the corresponding vibrational 
mechanisms. 
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FIG. 1: (Color online) Thermal equilibrium quantities for the configurational case. The dotted lines 
indicate the OD temperature calculated from free energy calculations(Figl2|). The order parameter 
rj is shown for two of the four sublattices of the LI2 structure, the order parameters not shown 
exhibit an identical behavior to that of sublattice (3. 
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FIG. 2: (Color online) Helmholtz free energy and entropy versus temperature for the configurational 
case. In (a) the solid line and the dashed line correspond to single realizations of the quasi-static 
heating and cooling processes, respectively. In (b) solid line depicts the entropy obtained from 
smoothing the —d{Fconf)/dT data, dashed lines correspond to the coexistence region, and the OD 
transition temperature T'^'^^ = 1925 ± 30K is estimated from the center of the coexistence region. 
The inset shows the entropy of a typical single realization where the coexistence behavior can be 
observed. 
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FIG. 3: (Color online) Behavior of the order parameters and entropy as functions of temperature 
in a logarithmic scale for the configurational case. In (a) the long-range order parameter rj vanishes 
at the OD transition, in contrast to the short-range order parameter a, which vanishes only at very 
high temperatures. This behavior is reflected in the configurational entropy showed in (b). The 
dashed line in (b) indicates the value of Sconfioo). 
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FIG. 4: (Color online) Behavior of the short-range order parameter o" as a function of temperature 
in heating and cooling numerical experiments at a rate of 0.02 K/MCS. From left to right, the 
three pairs of curves are for the fully relaxed, partially relaxed, and unrelaxed cases. The dashed 
lines indicate the OD transition temperatures, obtained from the crossing of the free energy curves, 
which are t/J = 1339 ± 20 K, TjJ = 1497 ± 40 K, and T^^ = 1635 ± 60 K. The curves are the 
smoothed data from averages over 10 samples. 
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FIG. 5: (Color online) Free energy and entropy as functions of temperature for the unrelaxed case. 
The dashed line indicates the OD temperature of T^J = 1635 it 60 K obtained from the crossing 
of the free energy curves. The lines in the entropy plot are obtained from the smoothing of the 
—d{Fur)/dT data. 
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FIG. 6: (Color online) Entropy as a function of temperature for the fully relaxed case (solid curve) 
and perfectly ordered LI2 phase (squares). The open circles depict the fully relaxed vibrational 
entropy, which is the difference between the fully relaxed total entropy and the fully relaxed config- 
urational entropy. The vertical short dashed line indicates the fully relaxed OD temperature. The 
open circles represent the smoothed data from averages over 10 samples. The solid and long dashed 
lines are fittings to the smoothed data and the error bars are smaller than the symbols (squares 
and circles). 
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FIG. 7: (Color online) Behavior of the total entropy as a function of temperature for the liquid 
phase and all cases studied. The inset shows the configurational case. The curves are the smoothed 
data from averages over 10 samples. 
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